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Abstract: We investigate with the help of analytical and numerical methods the reaction A + A — ► 
A on a one-dimensional lattice opened at one end and with an input of particles at the other end. 
We show that if the diffusion rates to the left and to the right are equal, for large x, the particle 
concentration c(x) behaves like A s x _1 (x measures the distance to the input end). If the diffusion 
rate in the direction pointing away from the source is larger than the one corresponding to the 
opposite direction the particle concentration behaves like ^4 a x -1 / 2 . The constants A s and A a are 
independent of the input and the two coagulation rates. The universality of A a comes as a surprise 
since in the asymmetric case the system has a massive spectrum. 
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1 Introduction 



The study of one-dimensional reaction-diffusion models far from thermal equilibrium is a field of 
growing interest. The dynamics of these models is characterized by non-trivial correlations so that 
ordinary mean field techniques fail. Therefore theoretical descriptions have to take local fluctuations 
into account. In general this is a very difficult task and approximation techniques are needed. 
However, there is a small number of exactly solvable models where we can derive exact results. 
The great interest in solvable models comes from the fact that their physical properties appear 
also in many other more complicated models. If we consider one-dimensional two-state models (i.e. 
models with only one species of particles), there are only two classes of exactly solvable systems. 
The first class includes diffusion (or exclusion) models which are solvable by means of Bethe ansatz 
techniques. The second class contains mainly models with an underlying theory of free fermions. 
The most important representative of this class is the so-called coagulation model which is the 
subject of the present work. 

Coagulation models describe particles which diffuse stochastically in a ci-dimensional space. 
When two particles meet at the same place they coalesce to a single one (A + A — * A). For a 
possible experimental realization see [Q. The theoretical study of coagulation models has a long 
history. It started with the observation that the critical dimension of the corresponding field theory 
is d c = 2 ||]. A breakthrough towards the exact solution of the one-dimensional coagulation model 
on a lattice was the introduction of so-called interparticle distribution functions (IPDF's) For 
certain reaction-diffusion processes the IPDF formalism leads to a hierarchy of decoupled differential 
equations similar to those obtained for the Glauber model Q . The most general conditions for the 
decoupling to occur as well as the cases when the underlying Hamiltonian can be diagonalized 
in terms of free fermions are given in Ref. ||. A variety of exact solutions were found for the 
coagulation model with or without back reaction (decoagulation A —* A + A) |H|. 

The one-dimensional coagulation model with spatial homogeneous external particle input at all sites 
has been studied extensively M and algebraic relaxation times have been observed. In this paper we 
investigate the same model with open boundary conditions and localized particle input at the ends of 
the chain. This problem was considered for the first time in Ref. jjj where various scaling relations 
could be obtained. In the present work we compute the particle concentration in the stationary 
state. We give the full solution on the lattice and in the continuum. The main motivation of this 
paper stems from the observation that in the mean field approximation (to be reviewed later) the 
density of particles has an algebraic decay and one can look if one has universality properties. 

The coagulation model studied in this paper is defined as follows. Particles of one species diffuse 
stochastically on a linear one-dimensional lattice with L sites. The diffusion may be biased due 
to some external force. If two particles meet at the same site, they coalesce to a single one. In 
addition particles are added stochastically with a given probability at the endpoints of the lattice. 
We use random sequential updates, i.e. we assume continuous time evolution which is described by 
a linear master equation. Altogether the dynamics is defined by the six nearest neighbor processes 
with rates shown in Figure |]. 

The master equation for a lattice of L sites can be written in a compact form by (for notations c.f. 



where the vector \P(t)) denotes the probability distribution. H is the time evolution operator which 
can be written as a sum of nearest-neighbor reaction matrices H nn +\ plus two further matrices for 



Ref. |) 




H\P(t)) 




1 



-e- 



• e— 



-• e- 



-Q- 



\ an 



biased diffusion 



-e- 



-• •- 



biased coagulation 



o- 



PL 



\ PR 



particle input at the boundaries 



Figure 1: Bulk and boundary processes in the coagulation-diffusion model 
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In the canonical basis of particle configurations ( |00), 1 0^4) , |^40), | ^4^4) ) these matrices read 



Hn.n 



+1 



( 

a L 


\ 



\ 

-an -c R 

aL an -cl 

cl + cr I M 

1 / n,n+l 



(1.3) 



and 



II 



PL 
-PL 



1 

It is useful to introduce some notations 
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c(i) denotes the particle concentration at the site i. When we will consider the continuum limit 
(the lattice spacing A — > 0), we will use the notation 
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Before giving our results we remind the reader of the improved mean field calculations in Ref.|7|. 
Assuming that the coagulation rate s is proportional to the concentration s = Qp{x), for large 
values of x (we take pr = and the source is at x = 0), the densities are 



q < 1 (bias to the left) : 
q = 1 (symmetric diffusion) 
q > 1 (bias to the right) : 
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Notice the algebraic fall-off for q > 1. 
The IPDF method is applicable if 
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These conditions are equivalent to r = t and s = 2. This will also be called the fermionic 
case (the Hamiltonian can be diagonalized in terms of free fermions ||, J| [K|). The special case 
q = 1 ,pr = pi = oo was already studied in a paper by Derrida et al [jll[ which was a source of 
inspiration for the present work. We list now our main results in the pr = case: 



a) Lattice in the thermodynamical limit (j fixed, L — > oo) 
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3 = 1 (1-9) 
9 > 1 (1.10) 
g<l (1.11) 



Several exact values for c(j) are given in Appendix A. 



b) Continuum and thermodynamical limit (x fixed, L — > oo) 
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c) Continuum and scaling limit [z 
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As a by-product, the one-hole functions are also obtained. This result is not trivial since they 
represent two-point correlation functions. 

In the long Section |2| and in the Appendices A and B these results are derived. Many of our 
calculations are extensions of results obtained in references Q and [ 10 1 from which we borrow the 
notations. Also in Section |2| we give the spectrum of the Hamiltonian in the one-hole sector. For 
q 1 the spectrum is massive in spite of the algebraic behavior seen in Eqs. ( 1.1C| ) and (|1 . 1 3|) . For 
q = 1 most of the excitations are massless (they coincide with those of the open chain (pi = 0)) 
but there are also some massive excitations with a mass given by pl. In many systems time-like 
and space-like properties seem to be coupled in the sense that long range correlations in time imply 
long range correlations in space. This is not necessary valid for stochastic models which are not 
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isotropic. As we will see in this model one can have short range correlations in time but long range 
correlations in space (for another example, see Ref. ||). 

In Section ^ we consider the problem of the universality: The coefficient - for the leading 
contribution in the q = 1 case (see Eqs. ( |Q| ) and ( 1.12j )), the coefficient y^pq^y^ for g > 1 (see 



Eqs. ( |1.10| ) and ( 1.13|) ) and the finite-size scaling function ( p.. 15 ). For this purpose we keep the 



definition y = <7 but leave the coagulation rates cr and cl arbitrary. For the open chain (no 
input) the spectra are known to be massless for q = 1 and massive for q > 1 |||. The modifications 
introduced by the boundary terms are supposed not to change radically the picture. Using Monte 
Carlo simulations (the details are given in Appendix C) we show that indeed for several values of 
cl and cr the expansion coefficients as well as the finite-size scaling functions are universal. 

The reader not interested in lengthy calculations can skip Section ^ and proceed directly to 
Section |3[ 

2 Exact solution in the clr = cr = q, cll = cl = g 1 case. 
2.1 Finite lattice calculations 

In this Section we give the full solution of the coagulation model with particle input at the bound- 
aries. We use the IPDF formalism || in which the whole problem is formulated in terms of 
probabilities for finding sequences of unoccupied sites (holes). In this basis the master equation 
leads to a hierarchy of sets of equations according the number of holes. It is known || that these 
sets decouple from the higher ones provided that the rates for diffusion and coagulation coincide 



(see Eq. fll.q) ). Therefore the one-hole sector decouples from the higher sectors and can be solved 



separately. In what follows we will assume that the above condition holds. 

For completeness we will consider the model with input at the left end (rate pi) and right end 
(rate pr). 4 We will actually show that by solving the problem with pr = one can obtain the 
general solution pa ^ 0. Although in principle feasible, we didn't look to the case when one has 
also output of particles at both end. 

Although not obviously needed for the study of the stationary state, we will also give the 
spectrum of the problem in the one-hole sector for two reasons. One is technical: the eigenfunctions 
and eigenvalues occur in the expression of the stationary state hole probabilities. The second one is 
related to the physical significance of our result: it is important to know whether one has massless 
or massive excitations. 

Let Q(j, m, t) denote the probability to find the sites j + l,j + 2, . . . m empty at time t. By a 
careful analysis of the elementary processes taking place at the edges of the hole one is led to the 
following equations of motion for the one-hole sector: 

• for holes which do not touch the boundaries (0 < j < m < L): 

^n(j,m,t) = qn(j-l,m,t) + q- l n(j + l,m,t) (2.16) 



4 In order to avoid confusion in terminology we give in parenthesis an alternative denomination used in the liter- 
ature: periodic boundary conditions (model on a ring), open boundary conditions (linear chain with closed ends), 
open boundary with particle input (chain with open ends). 
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+ qSl(j,m-l,t) + ^Ofem + M) - 2 (q + q l )Vt{j,m,t) 

• for holes touching the left boundary (0 = j < m < L): 

-^O(0, m,t) = q O(0, m-l,t) + q' 1 O(0, m+l,t) - (q + q' 1 + Pl) O(0, m,t) (2.17) 

• for holes touching the right boundary (0 < j < m = L): 

L, t) = qn(j - 1, L, t) + q-^j + 1, L, t) - (q + q' 1 + p R ) 0(j, L, t) (2.18) 

• for the hole extending over the whole chain (j = 0, m = L): 

-O(0, L,t) = -(p L +PR)n(0,L,t). (2.19) 

In these equations we have taken t) = 1. This leads to an inhomogeneous system of equations. 

Separating the time dependence and introducing rescaled probabilities 

n{j,m,t) = e- At q +j+m n(j,m) (2.20) 

we obtain the simplified system of equations 

faq + q- 1 ) - A)0(j,m) = 0(j-l,m)+0(j + l,m) + 0(j\m-l) + 0(j,m + l) (2.21) 

(q + q~ x +PL- A)o(0,m) = O(0, m - 1) + O(0, m + 1) (2.22) 

(q + q- l +PR-K)VL{3,L) = 0(j + 1, L) + Cl(j - 1, L) (2.23) 

(pL+PR-A)tt(0,L) = (2.24) 

with the inhomogeneous boundary condition 0(j, j) = g~ 2j . 

The homogeneous set of solutions describes the relaxational modes of the system. It is obtained 
by setting 0(j, j) = and can be computed easily by using similar techniques as in Ref. ||] which 
rely mainly on the invariance of the bulk equation ( 2.1 6| ) under reflections j <-> m and j «-> L — m. 
Denoting 

sinhf j arcsinhi(o + q' 1 — z)) 

9(3, z) = ) 2 - 4 ( 2 - 2 5) 

sinh ( L ar csinh ^ (q + q — 



the homogeneous solutions are 



$o(j,m) = g(L - j,PL)g(m,p R ) - g(L - m,p L )g(j,p R ) (2.26) 

(L), ftkj nkm 

$l'{j,m) = sm — g(L-m,p L )- sin— — g(L - j,p L ) (2.27) 

(R), nkj ixkra 

®k U, m ) = sm — g(m,p R )-sm—j—g(j,p R ) (2.28) 

nkj irlm irkm nlj 

®k,l{j,nT') = sin—— sm — — — sm — - — sin— — . (2.29) 
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They have the excitation energies 

A = Pl+Pr (2.30) 

A { k L) = q + q- 1 + PL- 2cos^ (2.31) 



L 

irk 

T 



= q + q' 1 + PR- 2cos^ (2.32) 



A M = 2(g + g" 1 ) - 2 cos^ - 2 cos^ (2.33) 

where 1 < k < I < L. In contrast to the coagulation model without particle input (pr = pi = 0) 
where the spectrum is massless for q = 1 and massive otherwise (q real), in the case of particle 
input, the spectrum is more complex. Even for q = 1 where most of the excitations are massless 
(Eq.( p.33]) ) we get some massive ones too (Eqs.( |2~3C| ) - ( |2.32| )). 

The derivation of the inhomogeneous (steady state) solution is more difficult. For symmetric 
coagulation on a ring with infinite particle input at a single site an exact solution has been found 



recently in Refs. [11, 12]. This solution applies to an open chain with symmetric diffusion (q=l) 



and infinite particle input at both ends (pi = pr = oo). It is given by j) = 1 and 
~, . 8 ^=!' sin z£ sin 4 (sin sin idp. _ sin e^r sin nU.\ 

n(j,m) = £ ^ h £ r ^ r-^ (j < m) 2.34 

L2 fe fc (cos£-cos^)(2-cosf-cos£) U ) K } 

where the prime indicates that the sum runs only over even values of k and odd values of I. Formally 
this solution can be expressed in terms of the two-particle excitations ( |2.29| ) by 

fittm) = ± £ Ai±M^!> (2 .35) 



k,l=i k ' 1 



where 



hi = \ (1 - (") fc+i ) ""J 8 " 11 ! ■ (2-36) 
2 cos ^ - cos ™ 

plays the role of a structure function. However, we one can prove that the general stationary 
solution for the asymmetric diffusion and finite particle input rates has the same structure and 
differs only in the structure function fkj. This function can be derived as follows. Let us symbolize 
a contraction of two functions over momentum indices k, I by {-,-)k,i an d similarly a contraction 
over spatial indices by (•, -)j, m - Then Eq. fl2.35| ) reads n = (/, x)fc,2 an d therefore the application 
of the discretized Laplacian A$ = A<3? yields AQ = (/, &)k,i which is zero everywhere except at the 
boundaries. Using the orthogonality relation <I>/v i>)j m ~ 5k k'$i V one can therefore compute 

/by 

/ ~ (/, 6S) ktl ~ (/, (*, $) j>m )k,i ~ ((/, ~ (AO, *> iim . (2.37) 

Carrying out these contractions it turns out that the structure function / consists of three parts 

fk,l = fi? + /£? + Af ■ (2-38) 

The first part /tj describes the asymmetric coagulation model with infinite particle input rates 
at both ends. It is given by 

Aoo) ( _ )k+ l a -2 L) (g + g- 1 ) 2 ginj sin ^ sin ^1 sin ^ 

^ " C U " V + ^-2cos^)(^ + ^-2cos^)' {2 - 69) 
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For q — > 1 (symmetric diffusion) this expression reduces to Eq. ( 2. 36] ), The other two parts depend 
on the input rates pl,Pr and read 



Jk,l 



AR) 
Jk,l 



■nk 



wl 



wl 



sin — sin -£ [cos 



cos^) 



2 (q + q- 1 + PL ~ 2 cos ?f) (q + q- 1 + p L - 2 cos f ) 



g-2L sin sin ^ ( cos *J 



cos 



2(g + ^ x +pr 



(2.40) 



(2.41) 



2 cos f) (q + q" 1 + p R - 2 cos f) 

The inhomogeneous solution of the difference equations ( |2.21| )-( |2.24 ) is then obtained by inserting 
Eq. (|2~38|) into Eq. ( ^35| ). 

For a fixed value of the lattice length L the holes probability function depends on three param- 
eters: pl, PR and q. By reversing the ends of the lattice one sees that: 



,m,p L ,PR,q) = £l(L -m + l,L- j + l,p R ,p L ,q u ; 



(2.42) 



Due to ( 2.38 ) the holes probability function obeys the following rule: 

,m,pL,PR,q) = ^(j ,m,p L ,0,q) + 0(L - m + 1 , L - j + 1 ,p R , , q' 1 ) 



Q(j , m , , , q) 
(2.43) 

Therefore it will be sufficient to study systems with particle input at only one boundary. Using 
( 2.43| ) one can relate physical quantities referring to systems with particle input at both ends with 
the ones computed for systems for which pi or p R is 0. As an example, the particle concentration 
at site j 

c(j) = 1 - (2.44) 

can be written as a sum ( |2.43| ): 

c(j,PL,PR,q) = c(j,p L ,0,q) + c(L - j + l,p R ,0,q' 1 ) - c(j,0,0,q) (2.45) 

Here c(j, 0, 0, q) is the particle concentration in the stationary state for input rates 0, i.e. it is the 
particle concentration of one random walker occupying the whole lattice ( |2.89| ) . 



2.2 The thermodynamic limit 

The formulas derived in the last Section are exact solutions for finite chains. We consider the 
thermodynamical limit. In this limit the right boundary is moved to infinity while the observer 
stays in a fixed distance j to the left boundary. We consider systems with no particle input at 
the right end {p R = 0). We are left with only two parameters, namely the input rate at the 
left boundary p = pl and the asymmetry parameter q. Carrying out the limit L — > oo in Eqs. 
(|2.38j )-( 2.41 ) one is led to a simple integral representation of the one-hole probabilities m). To 



this end it is convenient to introduce the quantities [i 



i(? + q- 1 - iz) ~\j\{q + q- 1 ~ iz) 2 - 1 (2.46) 



and its inverse 



/.J 1 = l( q + q -l-i z ) + yjl( q + r i- iz )2-l. (2.47) 

Using this notation, the one-hole probabilities in the thermodynamic limit are given by the elliptic 
integral 

^ ■») = 1 " i£ C iz (i - ?tf) - "* v -*) ' (2 ' 48) 
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A proof of this formula is given in Appendix A where also the particle concentration at the first few 
sites for infinite input rate is computed exactly. Let us now investigate the asymptotic behavior of 
the particle concentration 



c(j) 



1 



2iri q 



- z z 2 + p 2 ■ 



(2.49) 



for large j. Three cases have to be considered separately: 



i) Symmetric case (q = 1): 

For q = 1 the logarithm of the expression n z fi- z in Eq. ( 2.49| ) can be expanded to first order in z 
by 

log(/i*/i-*) = -yjM + 0(\z\V 2 ) (2.50) 

Rewriting the integral ( 2.49| ) by 

I r+co 



I r+oc j . 

ij) = ^— J dz exp[-y/2\z\ jjr(J,z) 



A/1 



r(j,z) = exp(y/2\z\ j) (jj, z fi- z ) 3 ( 



■ z z 2 + p 2 



) ("J 1 



(2.51) 
(2.52) 



and expanding r(j, z) in z the integral can be solved order by order. We obtain the series 



1 



C 0') = + T7 2 



1 



7TJ 7T/ 2 27TJ 4 \ 87T 7Tp 2 / j 



(2.53) 



This proves that in the fermionic case the first three terms in the large x expansion are independent 
of the input rate p. 

ii) Bias to the right (q > 1): 

In this case we find that the expression log(g 2 /jL z fj,— z ) in Eq. ( |2.49| ) can be expanded in first order 
by: 

logfa 2 ''fizfJ-z) = 



■ 2 1 + 1 z 2 - O(Z) . 



* (q 2 - 1)3 
Rewriting the integral ( |2.49| ) by 

c(j) = -^/* dze^(-q 2 ^±^jz 2 )s(j,z) 



s(j,z) = exp(g 2 J 2 + * 3 j z 2 ) (1 - ^^2) (gV^-*) 3 " 1 ~ M-l) 



2; z 2 + p 2 - 



(2.54) 

(2.55) 
(2.56) 



and expanding s(j, z) in 2 the integral can be solved again order by order. We obtain 

'3g 4 + 20g 2 -l (<? 2 -l) 3 \l 



1 



(<7 2 + 1)ttj 



1 + 



8(g 4 -l) 2q 2 (q 2 + l)p 2 I j 



-5/2n 



(2.57) 



We notice that, as opposed to the symmetric case, only the leading term is independent of the 
input rate. 

iii) Bias to the left q < 1: 

If the particles hop preferentially to the left, they accumulate at the left boundary and thus we 
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expect an exponential decay of the concentration profile. In fact, as can be seen from Eq. 



the expression q 2 i 1 c(j) is invariant under the replacement q — > q . This means that for a bias 
directed towards the left boundary the concentration profile decays like g 4j j -1 / 2 : 



^) = ^'U7^Tw7 + °(r 3/2 ) (2-58) 




We would like to remark that the series presented in this Section are asymptotic series since 
they are derived from an elliptic integral. 



2.3 The continuum limit 



An alternative way to describe the physics of the coagulation model with an external input source 
is to consider the continuum limit of the one-hole equations. This can be done by taking the lattice 
spacing A — > while keeping the two quantities 



2(q-q- 1 ) 



and 



P 



2 PL 



(2.59) 



constant. We then replace the empty-hole probabilities Q(x,y) by their continuous counterparts: 

(2.60) 

It is useful to rescale the hole density function taking Cl(x, y) = £l c {x, y)e~^ x+y ' > which verifies the 



n c (x,y) = n( J -,j) 



equation (see Eq, (|2.21|) ) 



A- -J n(x,y) 



(L > y > x > 0) 



(2.61) 



By solving the continuous counterparts of Equations ( 2,17| )- (|2.19 ) we determined the value of the 
holes density function on the boundaries. The solutions are: 



• Along the left boundary [x = , < y < L): 

sinh ( (L-y) 

^ if p ^ oo 

if p = oo 



sinh (l \J^+p) 





(2.62) 



Along the upper boundary (y = L , < x < L): 



Cl(x,L) 



' sinh 



x/L 



if f ^0 
if f =0 



On the diagonal (0 < x = y, < L ) the normalisation condition is: 

£l(x,x) = e~ rx . 



(2.63) 



(2.64) 
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Equation (|2.61| ) together with the boundary conditions ( 2.62| )- (|2.64| ) define a Dirichlet problem for 
the function y) . The formal solution is 



n(x,y) 



d 

dsn(x',y')—g(x,y,x',y' 
C on 



(2.65) 



where C is the contour along the boundaries, ^- the normal derivative and Q{x, y, x' , y') the Green 
function defined by 



(A - y)a(>, y, x', y') = 5(x - x')5(y - y') 

£(0, y, x',y') = Q{x, y, 0, y') = Q(x, L, x', y') = Q(x, y, x' , L) = 
Q(x, x, x', y) = G(x, y, x',x') = 



(2.66) 

(2.67) 
(2.68) 



The computation of the density function of the hole probabilities requires the computation of the 
Green function. This was done in two steps. First notice that: 



Q(x,y,x',y') = Q n (x,y,x',y') - G n (y,x,x',y' 



(2.69) 



where Gn(x, y, x' , y') is the Green function of the Dirichlet problem defined in the interior of the 
square < x < L , < y < L. Gn can be constructed easily by using reflection techniques. All 
what one needs to know is the Green function of the Dirichlet problem defined on the entire plane 
(with boundaries at infinity). We denote the last mentioned function with g(x, y, x' , y'). Summing 
up, we get: 



G(x,y,x',y' 



+oo 

E «p E 

a,/3=±l i,j=— oo 



g (^/(x - 2iL - ax'Y + (y - 2jL - /?y') 2 ) - (2.70) 



g(yj(y- 2iL - ax' f + (x - 2jL - f3y>) 2 



Once the density function Cl(x,y) is known one can compute expectation values of observables in 
the steady state. The local particle density is given by 



p(x) 



lim 



l-e^ x+ y^tt(x,y) 



y - x 



(2.71) 



y=x 



We also give a closed formula for the computation of the connected two-point function: 

G c (x, y) = (n x n y ) - (n x ) (n y ) (2.72) 
Here n x denotes the particle number operator at site x. Using the factorization properties of the 



two-holes probability function mentioned in [10, 11] it is easy to see that in the continuum limit 
we have: 



(2.73) 



2.3.1 The scaling limit in the symmetric case (q = 1) 

In the case of symmetric diffusion the differential equation (|2.61| ) reduces to a Laplace equation. 
The Green function can be obtained from ( 2.70| ) by replacing g(u) with ^ lnu. 
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Let us consider for simplicity first the case of an infinite particle input rate (p = oo). In this 
case the hole density function is zero for x = (see Eq. ( 2.62p ). We are left with: 



Io du Gix ' y ' x ' ,y,) l^=u~ L dx ' T^ g{x ' y ' x '' y ' 



y'=L 



Inserting Eq. ( |270| ) one is led to 

n» (*,!/) = 



(2.74) 



(2.75) 



+oo 



^ a (3 ^2 [arctan^ 

a,/3=±l ij=— oo 



a{2i - x/L) - (3(2j - y/L) 



(2i - x/Lf + (2j - y/Lf+a(2i - x/L)+(3(2j - y/L). 



arctan 



a + 2% — x/L 



+ 2j - y/L 



We see that the hole density function for p = oo depends only on x/L and y/L . From Eq. 
( 2.71]) we obtain the local particle density in the stationary state for infinite particle input rate: 



+oo 



(x/L-2i) + (x/L-2j) 



POO{X) vrL . J^ (x/L - 2if + (x/L - 2jf 



Defining z 

where 

*(*) = 



|j, Eq. (|2.76 ) can be rewritten as 

Lpoo(x) 



1 

7T 



E 



(z/2-i) + (z/2-j) 



(z/2-if + (z/2-jf 



(2.76) 



(2.77) 



(2.78) 



l ,j = — oo 

sinh(7rz) + sin(-7rz 



+oo 



cosh(7T2;) 



COS (TTZ 



sinh(2-7r(z/2 — k)) + sin(7rz) 



( 7rz ) \ cosh(27r(z/2 — k)) — cos(ttz) 



+ (k 



-k 



The function 3>(z), called scaling function, is odd and periodic with period 2. In the limit z — > 
the function diverges like (the dominant contribution is given by the first term in the second 
line of fl2.78|) ). For z — > 1, &(z) approaches the value 1, but the value itself in the point z = 1 is 
$(1) = . So the function is discontinuous for all integer arguments. 

We consider now the case of an arbitrary input rate p and look at the scaling regime (z fixed, L — > 
oo). If p is finite one picks up another contribution to the holes density function coming from the 
integration along the boundary segment (x = , < y < L) in ( 2.65 ). The difference between 
the values of the holes density function corresponding to infinite and finite input rates is: 



L dy' ni(0,y')^g(x,y,x',y') 



x'=0 



(2.79) 



We are here interested in the scaling and thermodynamical limit. For large values of the lattice 
length L , O|(0, y') behaves like exp(— y'y/p) ■ We expand the derivative of the Green function in 
Eq. ( 2.79 ) near y' = and get: 

f^(x,y)-^(x,y) 



16 +2? (x/L - 2i)(y/L - 2j) (y/L - 2jf - (x/L - 2if 

T7~ ^ 



7T 



I ,j = — OO 

1 sil^h(^/pL(l — u)) 
sinh (a/p L) 



(y/L - 2j)2 + (x/L - 2%Y 



u 6 + 0(u<) )du 
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It is easy to see that in the finite-size scaling limit one gets: 



**(z) = *(*) + Ti <S> con '(z)+0{- R ). 



(2.80) 



The first finite-size scaling correction function is 



+00 



$ corr (z) 



E (*/2 



(z/2 _ i)4 _ io(z/ 2 - i) 2 (z/2 - j) 2 + 5(z/2 - j) 4 



h3- 



(z/2-i) 2 + (z/2-j) 2 



(2.81) 



So we proved that in the scaling limit the particle density function scales. The scaling function 
<J?(z) given by equation ( 2.7£| ) is independent on the input rate. We get finite-size corrections of 
order L -4 for finite input rates. 



In Figure we show the function (solid curve) 



7T 7 

F(z) = — • *(*) 



(2.82) 



(with this definition F(0) = 1) together with finite-lattice calculations (see Eqs. (2.44) and ( |2.35| )) 
obtained for L = 2000 (p = 1) and for L = 800 {p = 00). We have given this figure for two reasons. 
First we observe that the scaling function has a nontrivial behavior at the opened end of the system 
{z oc 1). Next, although the scaling function was computed in the continuum limit, it applies to 
the lattice too. 



Let us consider the thermodynamical limit L — ► 00, x fixed. Using equations ( |2.71| ) and ( |2.70| ) 
(from the multiple sum of the latter one we take only the terms corresponding to i = j = 0) one 
can show that 



p(x) 



2 

TTX 



12 



Trp 2 x 5 



5040 

TTp 4 X 9 



+ o{p- b x - li ). 



(2.83) 



The asymptotic behavior (for large values of x) of the one-point function in the thermodynamical 
limit fl2~83p can be also obtained from the scaling behavior ( |2.80D , in the limit z — > . The result 
(2.82) is consistent with the expansion for the stationary concentration profile on a discrete lattice 
( 2.55 ). The only difference is that in the continuum limit all contributions with i > 2j + 1 
scale like A J ~ 2j ~ 1 and therefore vanish in the limit of vanishing lattice spacing A — > 0. 

We mention one last result concerning the connected two-point function. In the thermody- 
namical limit one can see (c.f. Ref. [O]) that the hole probability density for infinite input rate 
is: 



arctan 



7T 



(2.84) 



Using ( 2.73|) one gets the following expression for the connected two-point-function 



G c (x,x + d) 



16 1 



1 + v — v{2 + v) arctan [jt^ 



it 2 x 2 



(2 + 2v + v 2 ) 2 



(2.85) 



in the infinite input rate case. On the right hand side of ( 2.85| ) we denoted ^ with v. Notice again 
the algebraic fall-off. 
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2.3.2 Bias to the right (r > 0) 



Writing the Fourier transform of ( 2.66 ) one gets 



00 k-dk 



2- 



e- ikucose d9. 



(2.86) 



(2vr) 2 J k 2 + f 2 /2 

Thus the Green function of the problem defined on the entire plane is (see formulae (9.6.16) in [13 
and (6.532.4) in |L§]) 



9{u) 



2-k 



Kr 



-U 



(2.87) 



We use the standard notation Kj 



0, 1, 2... for the modified Bessel functions. 



In Appendix B we will prove that in the thermodynamical limit the asymptotic behavior of the 
particle density is given by 



p(x) 



2irx 



+ 



1 



4V27rf 



11 
~2 



1 



,3/2 



+ 0(x~ 5 / 2 ) 



(2.88) 



The leading term is the continuum correspondent of the one appearing in Eq. ( |2.57| ) and is inde- 
pendent of the input rate. 



We conclude this Section with some remarks concerning the influence of the right boundary on 
the particle density. 

As opposed to the symmetric case, one can check that for biased diffusion to the right the be- 
havior of the one point function in the thermodynamical limit is identical with the one in the scaling 
limit. One can give a qualitative explanation. In the stationary state, near the right boundary, the 
density function can be approximated with the one given by a single particle occupying the whole 
lattice (a random walker) 

Q- 2y TEf^r for q + 1 
c(y) = { ' (2.89) 



I 

L 



for q = 1 



where y = L — x. 



For biased diffusion to the right the density decays exponentially in y. So the influence of the 
right boundary is of short range and is not seen in the scaling limit. 

For q = 1 (or r = 0) the influence of the right boundary is much stronger. The density 
determined by a random walker near the right end of the lattice is constant (-^) and greater than 
the one obtained through the extrapolation of Eq. (2.53) (-3). This explains the linear behavior 



of the reduced scaling function for z — > 1 (see Figure 

One can also notice from the y and L dependence of c(y) that one has a characteristic length 
scale A = (21ng) _1 which is related to the inverse mass seen in the spectrum of the Hamiltonian 
(see Eq. ( |2.33j )). This observation is very interesting since it clarifies a puzzle which goes through 
this paper: how can a system with massive excitations in the time direction show an algebraic 
and, as we shall see in the next Section, universal behavior? The answer is that one looks at the 
concentration in the "wrong" way following the x dependence (away from the source) and not the 
y dependence (away from the open end). The discovery that this "wrong" way exists is probably 
the main achievement of this paper. 
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3 Numerical verification of the universality hypothesis 



In this Section we present the results of Monte Carlo simulations. The details are given in Ap- 
pendix C. We restrict the study to cases for which pr = and r > 0. We start with the case of 
symmetric diffusion (q = 1 or r = 0). First we look at the density profile (1 <C x <C L). As 
suggested by Eq.( |2.53| ), we fit the data by the function 



c(x) 



Ki 

x 



+ 



K2 



+ 



K 3 



x-' 



(3.1) 



When making the fits (by using the x 2 method) we took points x 6 [L/10 , L/2] in order to 
avoid finite-size effects. The estimates for K\ and K% for various input and bulk rates are given in 
Table |. The data presented here were obtained taking lattices of size L = 1000. 



s 


t 


V 


Ki 


K 2 


0.50 





1 


0.635 ± 0.001 


10.2 ±0.2 


0.50 





00 


0.639 ± 0.004 


9.8 ±0.8 


0.50 


-0.2 


1 


0.640 ± 0.010 


16 ±5 


0.50 


-0.2 


00 


0.640 ± 0.010 


13 ±3 


0.25 





1 


0.637 ± 0.002 


22 ±4 


0.25 





00 


0.632 ± 0.002 


24.4 ±0.4 


0.40 


-0.2 


1 


0.638 ± 0.002 


11 ±2 


0.40 


-0.2 


00 


0.642 ± 0.004 


10.4 ±0.6 



Table I: Estimates of the coefficients K\ and K% of the expression (3.1) for various input and bulk 
rates 

We notice that K\ is everywhere close to the value obtained in the fermionic case [s = 2, t = 0) 
namely K\ = 2/tt ~ 0.637. The values of K% are different if the bulk rates are different but as in 
the fermionic case they do not depend on the input rate. 

Since the leading term of the density profile is compatible with universality, one can go one 
step further and check if the scaling function <&(z) given by Eq. (|2.78j) in the fermionic case, is also 
universal. This function was obtained taking x/L = z fixed (x and L large): 



lim Lc(z,L) = $(z) 



(3.2) 



We define the function 



K(z,L) 



L c(z, L) 



1 



(3.3) 



which measures the deviation from universality and the finite-size effects. In Fig. || we give the 
data in the case p = 1 , s = 0.5 , t = for three lattice sizes. One notices that with increasing 
lattice size K(z, L) decreases, as it should. One should mention that for z = one expects K(z, L) 
to go to zero in the limit L — > 00 because of the universality of K\ in Eq.(|3.1|) and the uniform 
convergence of &(z) in the fermionic case. Thus the relative large values of K(z,L) observable in 
Fig. ^ for small values of z should not be a subject of concern. We have also done other simulations 
(not shown in Figure ^) for other input rates which show the same properties. 

In Figure ||, K(z, L) is shown for various input and coagulation rates for a lattice of length 
L = 1000. As one can see K (z, L) is small everywhere (for small values of z the convergence is slow 
but as mentioned above, for z = universality was checked already). 
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We now consider the asymmetric diffusion case (q > 1). As suggested by Eq. ( 2.57|) we fit the 
Monte Carlo data by the function 



c(x) = K[ 



/2-> 



-1/2 



3/2" 



-3/2 



(3.4) 



We choose q = y/2.5 (which corresponds to r = 0.857). 
K', ln = a/X ~ 0.369. 



in which case we get from Eq. ( |2,57j ) 
-l 



^ / 2 — y — 0-369. Notice that we have allowed a term ~ x 1 not present in Eq. ( 2.57| ). The 
data were collected for L = 1000 and the results are shown in Table pL 



s 


t 


P 




K[ 


1.000 


0.429 


1 


0.370 ± 0.003 


0.10 ±0.01 


1.000 


0.429 


OO 


0.368 ± 0.002 


0.11 ±0.01 


0.500 


0.214 


1 


0.368 ± 0.003 


0.90 ± 0.04 


0.500 


0.214 


oo 


0.369 ± 0.004 


0.82 ± 0.07 


0.250 


0.107 


1 


0.369 ± 0.002 


2.12 ±0.02 


0.250 


0.107 


oo 


0.369 ± 0.002 


2.09 ± 0.02 


0.361 


-0.181 


1 


0.368 ± 0.001 


1.27 ±0.02 


0.361 


-0.181 


oo 


0.369 ± 0.003 


1.25 ±0.04 



Table II: The coefficient Ky 2 and K[ of the expansion ( |3.4D for various input and coagulation rates. 
All data are for q = ^2l> (r = 0.857). 



As can be seen from this table the coefficient Ky^ is unchanged (universal). The K[ is in- 
dependent on the input, but depends on the bulk rates (similar to the Ki coefficient in Eq. 
(3.1)). Finally we notice that the values of K[ get smaller if we approach the fermionic case 
(r = t = 0.857, s = 2). 

To sum up, in the symmetric diffusion case the Monte Carlo data suggest that the large x 
behavior and the scaling function are universal: they are independent of the cl, cr and p rates. In 
the asymmetric diffusion case, the large x behavior is also universal. 



4 Connection with other models 

It is a well-known fact that the coagulation model A+A^A and the annihilation model ^4±^4 — ► 
belong to the same universality class. This equivalence is due to the existence of a local similarity 



transformation between their time evolution operators [15]. 



We now use this transformation in order to apply the results of the preceding Sections to a 
coagulation-annihilation model (called CA) with boundary effects which is defined by the following 
processes and rates: 

^40 — ► ®A diffusion to the right at rate clr 
$A — > A$ diffusion to the left at rate cll 
AA — ► ®A coagulation to the right at rate cr 
AA — » A$ coagulation to the left at rate cl 
AA — » 00 pair annihilation at rate k 

In addition particles are absorbed (desorbed) at rate 7 (5) at the left boundary. In the configuration 
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basis the time evolution operator H CA = + J2 n =l ^nn+i * s gi ven by 



h: 



CA 
n,n+l 



( o 








a-L -an 
-at a R 



-K 
-OLL 



\ 











t ca 



7 

-7 



(4.1) 



K + CR + CL J 



As shown in Ref. |T5[, the coagulation model (1.3) and the generalized annihilation model (4.1) 
are related by a local similarity transformation H CA = \J H coag XJ~ x 



U 



u (g) u <g> . . . (g) u 



II 



1 I- a 
a 



(4.2) 



where a is some parameter. The rates of the coagulation-annihilation model are related to those of 
the original coagulation model by 



dL,R 

cl,r 



7 
6 



a L,R 

CL,R + 
1 - a 



■( a R,< 



a L,R ~ Crl) 



a 
ap 

(1 - a)p 



(cl + cr) 



(4.3) 



Notice that if the original model had only input of particles the equivalent coagulation-annihilation 
model has both input and output of particles. 

Because of the simplicity of the transformation the n-point density-density correlation functions in 
the coagulation and coagulation-annihilation model are related by 



( T h T : 



\GA 



31 



n I \coag 
a V]\ T 32 ■ ■ ■ T 3nl 



(4.4) 



5 Conclusions 



In the present paper we investigated the coagulation-diffusion model with particle input at one 
boundary using both analytical and numerical methods. The results show that spatial long-range 
correlations play an essential role and that some physical properties are universal with respect to 
the input and the coagulation rates. 

We started our analysis with a simple space-dependent mean field approximation. It predicts 
algebraic behavior of the particle density in the stationary state for both symmetric and biased 
diffusion. However, rigorous results require an exact solution of the problem. To this end we solved 
the full problem by using the IPDF formalism. This formalism can be used only if the coagulation 
rates coincide with the diffusion rates, which corresponds to free fermions in the Hamiltonian 
language. The large x behavior (x is the distance to the source) of the particle density was computed 
in the thermodynamical limit both for the lattice and the continuum version and the results are 
compared. These painful calculations were done for symmetric and asymmetric diffusion. In the 
case of symmetric diffusion the scaling limit (x/L fixed, L is the lattice length) was obtained. 

Monte Carlo simulations show that the coefficient of the leading terms of the asymptotic ex- 
pansion of the density in the thermodynamical limit are universal: they are independent of the 
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input rates (this was to be expected from mean-field) and on the coagulation rates. The scaling 
function is also universal in the symmetric case. It is trivial in the asymmetric case (it coincides 
with the leading term of the large x behavior of the density). These results were to be expected 
from common sense in the symmetric case but not for the asymmetric case. The reason is the 
following one: the relaxation spectrum of the system is massless in the first but massive in the 
second case. There exists a myth according to which if there are lengths in the time evolution 
there should be lengths in the space correlations. A counter-example can be found however in the 
kinetic Ising model (see Ref. Q). In the coagulation-diffusion model the picture is more perverse : 
if one looks at the concentration starting at the opened end, one finds an exponential fall-off but 
an algebraic and universal behavior if we start at the source end. 

The message of this paper can be extended to the problem in which we add pair-annihilation 
in the bulk and an output of particles at the source. What is still missing is a proof of universality 
which goes beyond numerical checks. This can be done using field theoretical methods a la Cardy 
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A Proof of the solution in the thermodynamic limit 



In this Appendix we prove the integral representation for the one-hole probabilities in the 
thermodynamic limit (2.48): 



1 



2m 



J-oo V -Z 



z 2 + p 2 



(A.l) 



Instead of deriving this formula from the finite-size solutions (|2T3S|)- (|2.41D by taking L — ► co it is 
much simpler to prove that Eq. (pO|) is a solution of the one- hole equations Q2.16[ )-( p.l9|) . We first 
notice that il(x, x) = 1 because of the antisymmetry of the integrand. In order to verify the bulk 



equation ( |2.16| ) let us introduce the notation g(x,y,z) = q x+y (fj,^fjt 



Using Eq. ([HI) 



one can show that 

qg{x ~ hV,z) + <rV 2; + 1,2/,-z) +qg{x,y - M) + q~ 1 g{x,y + l,z) = 2(q + q' 1 )g(x,y,z) (A. 2) 



This relation implies that Eq. ( |A.l| ) satisfies the bulk equation fl2.16|) . The last step is to verify 
the left boundary condition Eq. ( |2.17 ). Obviously it is equivalent to prove that 



Kv) 



q fi(0, y - 1) + q- l n(0, y + l)-{q + q- 1 + p) fi(0, y) 



(A.3) 



-p + 



p 2 q y 



I d i( 



2-7T J z ^z + ip z — ip 
is equal to zero for all y = 1, 2, . . . oo. For y = we get 



Ho) 



+oo 



dz 



1 



z 2 + p 2 



(A.4) 
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For y = 1 one has to solve the integrals 

,2 



h(l) 



-p + 



QP 
2vr 



— oo 



Z 1 + p 2 



+ °° ip(n z — fi- 

z 



OO 



2 Z 



(A.5) 



by standard integration techniques in the complex plane. It turns out that all contributions cancel 
except at z = in the second integral so that h(l) = 0. Using Eq. ( |2,46| ) it is now easy to derive 
the recurrence relation 

q^hiy) +qh(y-2) = (q + q' 1 ) h(y - 1) (A.6) 
so that h(y) = for y = 2, 3, . . . oo follows by induction. This completes the proof of Eq. ( |A.l| ). 



Let us finally consider the case of infinite input rate where Eq. (A..1) reduces to 



£l(x,y) 



1 



qX+y 
2ni 



r+°° dz / 

J-oo Z V 



(A.7) 



This expression turns out to be a combination of elliptic integrals. This allows us to compute 
the particle concentration exactly although the expressions become very complicated as x and y 
increase. For example the particle concentration in the steady state at the first four sites is given 
by 



c(l) 
c(2) 

c(3) 
c(4) 



1 

2 

3n 
2 

157T 



(q 2 + 1)(1 + 6q 2 + q 4 )B - (q 2 + l)(q 2 - 1) 2 K 



2q 4 



(q 2 + 1)(4 - 15q 2 - 34q A - 15q 6 + 4g 8 )E 



(q 2 + l)(q 2 - 1) 2 (4 - 15q 2 + 4g 4 )K 



+ 3g 4 - 2q b 



1057T 



(q 2 + 1)(12 - 28q 2 + 45q 4 + 238q 6 + 45q s 
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+ 12g 12 )E 



(A.8) 
(A.9) 

(A.10) 
(A.ll) 



(? + 1)0? - 1) (12 - 28g 2 + 69g 4 - 28g b + 12<f )K - 3q b - Aq 



where 



E = 
K = 

are elliptic integrals of the first kind. 



tt/2 



tt/2 



d9 (] 



(q + q- 1 ) 2 
4 



(q + q- 



sm 



sm 



1/2 
-1/2 



(A.12) 
(A.13) 



B The one point function for biased diffusion to the right in the 
continuum limit 



In this Appendix we give a proof of equation Q2.88 ). We concentrate on the thermodynamical 



limit of models in which the particle motion is subject to a drift pointing away from the source 
which is situated at the left boundary (x = 0) and we are interested in the large x behavior of the 
density. The starting point is the contour integral ( p. 65 ). In the thermodynamical limit the Green 



function of the Dirichlet problem is defined by the i = j = term of the multiple sum of (ETrcD . 
We are left with two contributions to the holes density function 

n c (x, y ) = nSo(x.y) + n;(x, y ). (B.l) 
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The first one comes from the integration along the diagonal boundary half-line (0 < x' = y' < 

oo) 

nSo(*,v) = ^ E (*y - M T ^^e-^+yVV da (B.2) 



a,/3=±l 



where n = \J (a — 9 x +l 3 V ^j + ^ ° tx . This is the holes density function in the case of an infinite 
input rate. 

The second contribution comes from the integration along the left boundary half-line (0 < 
y' < oo , x' = 0) and is p dependent 

Ifl r°°( ^if% N ) — 

fi c (x } = y ae r(*+«)/2 / J y x v^; _ x v^; i e -\f¥ r i±H> a da (B3) 



where r2 = \J ( a ~ ax ) + U 2 and r 3 = J (a — ay^J + x 2 . 

The density profile is determined by the holes density function in the limit y — ► x (see Eq. 
( |2,71| )). The behavior of the integrands appearing in ( |B.2| ) and ( p.3[ ) is given by terms of the form 
K\(v)/u. For u —* the modified Bessel functions diverge like K v (u) ~ tt~ v (for Re(z/) > 0). 
The only dangerous term is the one containing r\ which vanishes for y — > x and a — > x when 
a = /? = 1. The corresponding term in flB.2|) 



2 



n§(x,y) = l^- g ) .r V/ v L e -na-(* + v)m da (B . 4) 



2 



27T /7 ^,\2 



determines the asymptotic behavior of the particle concentration, in the thermodynamical limit. 
We start with this term. We use the fact that K\ (y/u) / \/u is the Laplace transform of exp(— 1/ (4i)) 



(see Eq. (29.3.122) in ref. [|T3|]). The equation (B.4) can be rewritten 

nSfoiO = r 2V -^- ■ | o °°exp^-tf 2 (^)^ dt |^exp^-t(fG + l)^ da (B.5) 
After integrating over the variable a we get 

ns(.,») - i - ^ -r exp (- (2f2 ^) 2 ) erfc ^ " s) dt (B ' 6) 



Here erfc stands for the complementary error function. From ( 2.71 ) we get that the contribution 
of Qq(x, y) to the particle density is 

r f°° ( 1 \ 

p0 ^ = ' Jo eric ( tfx ~ Yt) dt ^ B ' 7 ^ 

With the change of variable 

t = " + ^ + 2 " X (B.8) 
2rx 



5 We note that the derivatives of the Bessel functions with respect to their argument are K (u) = —Ki(u) and 
K[{u) = -\{K {u) + K 2 {u)) 
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we get 



Po(x) 



1 



1 



\fuj 2 + 2fx 



Vu> 2 + 2fx 



erfc(a;) duo 



After integrating by parts one gets 

1 roo _ 2 

/)o(a;) = — / v w 2 + 2f s e w duo 
ttx Jo 

Expanding in powers of lo 2 /x we finally obtain: 



(B.9) 



(B.10) 



+ + 0(x" 5 / 2 ) 



2vrx 8\/27r? x 3 / 2 



Po{x) = 

The rest of the holes density function 

gives also a contribution to the particle density 



Prest\X) 



d_ 

dy 



(B.ll) 



(B.12) 



(B.13) 



In all the integrals contributing to ^.est-> * ne arguments of the Bessel functions (rj, i = 1,2, 3) are 
greater then x (or y). We can thus use the asymptotic expansions of K v {u) and replace these 
functions with exp(— u) ■ \Jtt/{2u) . One obtains an integral expression of p r est{x). Due to the 
exponential falloff of the integrands appearing in (B.13), one can expand them in powers of a/x 
(where a is the integration variable). After some computations one gets that for x — > oo p res t(x) 
decays like x -3 / 2 . 

Summing up we get the following formula for the asymptotic behavior of the density profile 



p(x) 



r 1 

+ 



2-7TX 



4V27rf 



11 

T 



P 



l 



,3/2 



+ 0(x 



-5/2n 



(B.14) 



We notice that the leading term is independent of the input rate p but the next to leading term is 
p dependent. 



C Details on the Monte Carlo simulations 

In this Appendix we explain how we have done the simulations of the coagulation-diffusion 
model with particle input at one boundary (pr = 0). We can simplify the notation and use the 
symbol p instead of pi for the input rate. 

The simplest way to simulate reaction-diffusion models is to use a Monte Carlo algorithm with 
random sequential updates. However, this algorithm is not very efficient for the present problem 
since particle densities are very low and therefore most of the updates take place at empty sites. 
This is why we used a different method in which the positions of the particles are stored rather 



than the occupation numbers of the sites (for details see [18] and references therein). This 'direct' 
method is much faster than the first one. It is defined as follows. At the beginning the lattice is 
empty. As long as the total number of particles in the system is the following steps are repeated: 
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• choose a randomly a number a, between and 1 

• occupy site 1 with a particle if a < pAt 

• leave the site unoccupied if a > pAt 

• increment the time t — > t + At 



The parameter At is the time discretization (see Ref. 
system, the 'direct' Monte-Carlo algorithm is started: 

• choose a particle at random and one of its neighboring sites 

• update the configuration of the chosen pair with help of a random number and by considering 
the bulk reaction rates 

• increment the time t — > t + ^ where N is the current total number of particles in the system 

• if site 1 is empty try to occupy it by comparing a random number with p^. If p = oo one 
keeps the site 1 occupied at all time-steps. 

In order to test the accuracy of the "direct" Monte Carlo method, we simulated some systems 
for which analytical data is available. The agreement is very good. For q = 1 we used a lattice 
of length L = 200 and took p = 1 and p = 0.01. Only for the first 10 — 15 sites the two sets of 
values for the density profile are slightly different (the relative difference is of less than 10% ). For 
the other sites the difference between the two measurements of the density profile is zero within 
numerical errors. We also compared data obtained for systems characterized by q = y/2.5, L = 20 
sites and p = 1 and p = 0.1. Although the lattice length is small, the two sets of measurements of 
the density profiles coincide for all sites within numerical errors. 

The quality of the Monte Carlo simulations is higher in the case where the particle diffusion 
is biased to the right in comparison with the clr < case. This has two reasons. On the one 
hand in the an > ai case the total concentration of particles in the stationary state is larger as 
compared to the symmetric case. It is more likely to reproduce through simulations a distribution 
with a higher total number of particles. On the other hand the relaxation of these system to the 
stationary state occurs much faster than in the clr = aL case since the time operator has massive 
excitations for an i= ol • Therefore less CPU time per run is necessary and thus the numerical 
errors of the measurements are smaller. 

It is a well known fact that the quality of the Monte Carlo determinations is limited by the 
accuracy of the random number generator. If the number of steps requiring random numbers is too 
high, at some point the generator produces correlated numbers. This limitation poses some prob- 
lems in the case of symmetric diffusion, for large lattices. This explains the unphysical oscillations 
of the data corresponding to L = 1000 in Figures || and ^ The choice of a better random number 
generator implies the increase of the CPU time needed to perform the simulations. 

Since we are interested in the stationary properties of the system we stopped each simulation 
run at a value of t = t max such that at least in the time interval [t max /3 , t max ] the average total 
number of particles is fluctuating around a constant value. We used a double averaging technique. 
We took 100 equidistant time points between [0.9 • t max , t max ] and measured our observables in 
each of them. For each Monte Carlo run of the program we got a preliminary value by averaging 



18]). After the first particle entered the 



21 



over this 100 determinations. Afterwards we averaged these preliminary values over all MC runs. 
The number of runs performed for each system was between 4 and 50 thousand, depending on the 
lattice length. Due to CPU time limitations, the number of runs performed decreases with the 
lattice length. 

For the data presented in Figs. [| and |] we used coarse-graining for obvious reasons, this is 
reflected in the horizontal error bars. 
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